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1 .  INTRODUCTION 

The  maximum  entropy  method  (MEM)  has  heen  successfully  used  (ref , 1 ,3)  to  estimate 
the  power  spectrum  of  severely-truncated  time  series.  The  enhanced  resolution 
of  the  MEM  power  spectral  estimates  are  due  to  the  estimation  of  a  set  of  linear 
prediction  filter  coefficients  which  can  be  used  to  extrapolate  either  the  data 
or  the  covariance  function  of  the  data  outside  a  finite  observation  interval. 

In  frequency  domain  beamforming  the  estimation  of  the  wavenumber  power  spectrum 
using  data  from  a  linear  array  of  equispaced  sensors  is  well-known  to  be  analogous 
to  the  estimation  of  the  power  spectrum  of  a  time  series.  As  shown  in  Section  2 
the  wavenumber  power  spectrum  at  each  frequency  of  interest  is  the  Fourier  transform 
of  a  complex  spatial  covariance  function  which  can  readily  be  derived  from  the 
cross-power  spectral  matrix.  For  an  array  of  limited  aperture  MEM  can  then  be 
used  to  extrapolate  this  spatial  covariance  function  independently  at  each  frequency. 
It  should  be  pointed  out  that  the  estimates  of  the  frequency  wavenumber  spectrum 
discussed  in  this  report  differ  from  the  two-dimensional  MEM  recently  proposed  by 
Wo6ds(ref .4) .  In  the  application  considered  here  the  extrapolation  in  the  time 
domain  is  unnecessary  since  increased  frequency  resolution  can  easily  be  obtained 
by  increasing  the  length  of  the  time  data  sequence  from  each  sensor. 

In  Section  4  the  estimation  of  the  spatial  covariance  function  of  data  from  a 
linear  array  of  equispaced  sensors  is  discussed  and  the  effect  on  the  MEM  of  biases 
in  estimates  of  the  covariance  function  due  to  windowing  is  shown. 

Finally  in  Section  5  the  methods  discussed  in  the  previous  sections  are  used  to 
estimate  the  maximum  entropy  frequency  wavenumber  power  spectrum  of  sonar  data 
from  a  linear  array  of  equispaced  sensors.  For  comparison  the  Fourier  estimates 
of  the  frequency  wavenumber  power  spectrum  were  calculated  and  are  also  presented. 

The  sensitivity  of  the  MEM  to  the  length  of  the  estimated  spatial  covariance 
function  is  illustrated  by  means  of  an  example. 

This  work  is  part  of  a  continuing  R§D  programme  in  signal  processing  for  underwater 
acoustics  and  has  been  carried  out  under  task  DST  79/069. 


2.  THE  FREQUENCY  WAVENUMBER  POIVER  SPECTRUM 

In  beamforming,  the  power  spectrum  of  noise  incident  upon  the  array  is  often 
estimated  as  a  function  of  frequency  and  the  angular  co-ordinates.  Because  of 
the  rotational  symmetry  of  a  linear  array  only  the  total  power  density  incident 
on  the  array  from  a  cone  defined  by  an  angle  9  relative  to  the  axis  of  the  array 
can  be  measured.  Thus  the  noise  field  can  be  estimated  only  as  a  function  of 
frequency  and  the  angle  9.  An  alternative  approach  (discussed  in  more  detail  in 
reference  5)  is  to  use  the  single-dimension  co-ordinate  wavenumber*  k,  defined  by 

k  =  sin  9 /\ 

instead  of  9.  The  rationale  for  using  the  frequency  wavenumber  approach  has  been 
discussed  in  reference  5  and  is  adopted  in  this  report. 


sampling  interval,  it  follows  that  if  A(f,k)  is  defined  by 


*  sin  9 fK  is  strictly  k^,  the  component  of  the  wavevector  along  the  axis  of  array. 

Since  only  one  dimension  is  considered,  the  x  subscript  is  dropped  and  k  is 
referred  to  as  wavenumber. 
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M/2  N/2 


A(f,k) 


lim  lim  1  1  V”"  ■>  -TJtiflT 

Nt*oo  NHoo  M+1  N+1  y  ^  ^  ®  ° 


,  (13 


j=-M/2  <=~N/2 


where  d  is  the  separation  of  adjacent  receivers,  then  the  two-dimensional  power 
spectral  function  is  given  by 

S(f,k3  =  ^(f,k)  A*(f.k)  >  , 

where  O  denotes  ensemble  averaging. 

Defining 


x.(f3 


lim  1 
N»’«>  N+1 


.  -2Jrif£T 
x^(£T^3  e  o 


I  =-N/2 


M/2-j 


r^Cfj 


m  ^  <511  cf)  «.ii 


i=-M/2+j 


then  it  can  easily  be  shown  that 


S(f.k)  -  1^1  HTl 


j=-M/2 


By  analogy  with  time  series  analysis  it  is  natural  to  term  r^ (f )  the  spatial 

covariance  function  since  its  infinite  Fourier  transform  is  the  wavenumber 
spectrum.  One  important  difference  is  that  r ^ (f)  is  in  general  complex  whilst 

for  time  series  analysis  the  covariance  function  in  general  is  real. 

For  an  array  of  finite  extent  r^ (f)  can  only  be  estimated  for  I  jl  <  M'  where  M'd 

is  the  aperture  of  the  array.  Hence  a  natural  estimator  is  obtained  by  replacing 
the  infinite  sum  in  equation  (1)  with  a  finite  one.  For  an  array  with  a  finite 
number  of^receivers  and  for  a  finite  number  of  sampled  data  points  the  Fourier 
estimate  Sg(f,k3  is  defined  to  be 


Sg(f.k3  =  <  Ag(f.k)  A*(f,k)  > 
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where 


Ag(f.k) 


M/2 

1 

M+1 

j=-M/2 


X.(f)e 


-27rikjd 


and  where 


x.(f) 


J_ 

N+1 


,  -TkxUt 
x.Ut  )e  o 

J  0-’ 


^=-N/2 


(3) 


The  subscript  B  has  been  used  to  denote  that,  in  forming  this  two-dimensional 
periodogram,  a  Bartlett  window  has  been  implicitly  applied  to  the  spatial 

covariance  because  of  the  presence  of  the  term  in  the  summation. 


3.  MAXIMUM  ENTROPY  ESTIMATES  OF  THE  WAVENUMBER  SPECTRUM 

At  each  frequency  the  MEM  may  be  used  to  extend  the  spatial  covariance 
function  derived  from  a  finite  aperture  array.  As  discussed  in  the  introduction 
each  frequency  is  treated  separately  on  the  assumption  that  time  series  outputs 
from  all  receivers  can  be  made  sufficiently  long  to  achieve  any  desired  frequency 
resolution.  Thus  the  method  of  Woods  for  estimating  two-dimensional  maximum 
entropy  spectrum  need  not  be  used  and  a  simple  extension  of  the  method  of  Edwards 
and  Fitelson(ref . 6)  into  the  complex  domain  can  be  used.  The  method  is  outlined 
in  this  section  for  the  sake  of  completeness. 

The  entropy  at  any  given  frequency  f  for  a  process  with  a  wavenumber  spectral 
density  S(f,k)  may  be  defined  as 

H(f)  =  /j,  log  S(f,k)  dk  (4) 


where  the  integration  runs  over  the  values  of  k  for  which  S(f,k)  is  non  zero. 

The  maximum  entropy  estimates  §,,„(f,k)  of  SCf,k)  are  derived  by  maximising  (4) 

Mt  ^ 

subject  to  the  constraint  that  the  inverse  Fourier  transform  of  Sj^(f,k)  is  equal 
to  (f) .  (A  more  recent  derivation  of  the  method  by  Newman(ref , 7)  does  not 
require  equality  but  a  close  approximation  to  the  observed  r^Cfj's.  This  method 
is  not  considered  in  this  report). 

The  problem  of  maximizing  (4)  subject  to  the  constraint 

k 

o 

I  dk  =  ?.(f) 

-k 
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-  4  - 


for  j  =  0,  ±1,...,±(N-1) 
and  where 


f!(f)  = 

can  readily  be  solved  by  introducing  Lagrange  multipliers  a  la  Edwards  and 
Fitelson.  The  solution  is 


It  is  worth  commenting  that  Pj^  can  be  related  to  final  error  of  a  prediction 

filter  which  uses  N-1  previous  values  to  estimate  one  lag  ahead.  A  wealth  of 
papers  exists  on  efficient  algorithms  for  solving  the  Toeplitz  equation  (6) 
(see  reference  8) . 

As  a  final  comment  reference  should  be  made  to  an  alternative  method  developed 
by  Burg(ref.2)  whereby  the  prediction  filter  coefficients  (ie  the  ci’s)  can  be 
estimated  directly  from  the  data  without  the  need  to  estimate  the  covariance 
function. 


4.  ESTIMATION  OF  THE  SPATIAL  COVARIANCE  FUNCTION 

The  application  of  the  MEM  described  in  the  previous  section  requires  that  the 
spatial  covariance  function  at  each  frequency  of  interest  be  estimated.  From 
equation  (2)  a  reasonable  estimate  of  r.(f)  would  appear  to  be 


-1 
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M/2-j 

^  Vi'«  ^ 

i=-M/2+j 


where  <>„  denotes  averaging  over  P  independent  estimates,  ie 


<Xi(f) 


^  fki 

X.^.  (f) 


t  ^  m 


In  general  a  good  approximation  to  equation  (8)  can  be  obtained  by  segmenting 
the  total  time  series  into  P  consecutive  blocks. 

It  can  easily  be  shown  that  the  estimate  r^.  (f)  given  by  equation  (7)  is  biased. 

The  process  of  averaging  reduces  the  variance  of  the  estimator  r.(f)  but  not 

1  ^ 

the  bias  which  is  caused  by  the  implicit  weighting  function  (ie  the  Bartlett 

weighting  in  equation  (7)).  To  illustrate  the  seriousness  of  this  type  of  bias 
in  rj (f)  consider  the  example  of  sine  wave  in  white  noise.  The  time  series 

argument  has  been  chosen  here  but  it  can  easily  be  shown  that  it  is  trivial  to 
extend  the  following  argument  to  the  spatial  domain.  Replacing  the  finite 
average  by  the  ensemble  average,  equation  (7)  reduces  to 


Is.  +  a  cos  2wfjT 
M+1  JO  0 


where 


1  j  =  0 
0  j  ^  0 


and  o,  is  the  signal-to-noise  power  ratio.  A  plot  of  r^ (f)  is  given  in  figure  1 

where  the  Bartlett  (or  triangular)  weighting  over  the  interval  (0-T)  can  readily 
be  seen  for  a=lO.  The  Fourier  spectrum,  the  maximum  entropy  extrapolation  of 
the  covariance  function  and  the  corresponding  maximum  entropy  spectrum  are  also 
shown  in  figure  1. 

As  can  be  seen  the  MEM  interprets  the  weighting  on  the  covariance  function  as  a 
modulation,  extrapolates  both  the  modulation  and  the  carrier  and  consequently 
splits  the  single  spectral  line  into  two  closely  spaced  lines  beating  together. 
Thus  the  bias  introduced  by  using  equation  (7)  as  an  estimator  of  the  covariance 
function  has  a  radically  different  (and  vastly  misleading)  effect  on  the  MEM 
power  spectral  estimates  than  it  does  on  the  Fourier  ones  where  it  permits  a 
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well-known  trade-off  between  main  lobe  and  side  lobe  responses.  It  should  be 
emphasized  that  this  is  not  a  weakness  of  the  MEM  -  given  the  distorted  covariance 
matrix  it  made  optimum  use  of  the  available  information.  The  same  argument  is 
directly  applicable  to  the  spatial  domain  where  the  corresponding  effect  would  be 
the  splitting  of  a  single  arrival  into  two  distinct  ones. 

The  technique  adopted  in  this  report  for  removing  the  bias  is  to  replace 
equation  (7)  by 


M/2-j 


2M+l-j| 


<  X.  (f)X*  .(f)  >  _ 
1  ^  j  +  i  ^  P 


-M/2+j 


(9) 


which  is  now  an  unbiased  estimator  of  r.(f). 
estimator  of  r^ (f)  increases  linearly  with  I  jl 


Unfortunately  the  variance  of  this 
.  Thus  the  variance  of  the  estimate 


of  greatest  lag  will  be  M  times  the  variance  of  the  zero  lag  estimate.  An  effect 

of  statistical  variations  in  estimates  of  the  r^ (f)  is  move  the  position  of  poles 


in  the  response  function  of  the  linear  prediction  filter  coefficients.  In 
particular  the  effect  of  noise  on  poles  just  inside  the  unit  circle  can  be  to 
push  them  outside  the  unit  circle  and  consequently  the  method  becomes  unstable. 
This  results  in  the  maximum  entropy  extrapolation  expanding  instead  of  decaying 
exponentially.  As  discussed  in  the  following  section  this  can  be  avoided  by  only 
using  M'  (where  M'  the  available  lags. 

4 


5.  APPLICATION  TO  REAL  DATA 

The  method  described  in  the  previous  sections  has  been  applied  to  the  estimation 
of  the  frequency  wavenumber  power  spectrum  of  acoustic  data  from  a  linear  array 
of  32  equispaced  hydrophones. 

A  block  diagram  of  the  processing  is  shown  in  figure  2.  The  outputs  from  each 
hydrophone  were  narrowband  filtered  (by  means  of  the  EFT  routine)  into  40  indepen¬ 
dent  frequency  bins.  The  algorithm  is  more  efficiently  implemented  by  inter¬ 
changing  the  order  of  summation  in  equation  (9),  ie  the  cross-power  spectral 
matrix 


X.(f)  X*^.(f) 

was  estimated  at  each  integration  and  then  was  averaged  along  the  diagonals  to 
give  an  estimate  (ie  equation  (9))  of  the  spatial  covariance  function.  This 
estimate  was  then  smoothed  by  averaging  over  300  independent  samples  of  the 
spatial  covariance  function.  The  assumption  was  made  that  the  frequency  wave- 
number  power  spectrum  remained  stationary  over  this  period. 

At  each  frequency  lags  of  the  spatial  covariance  function  were  used  to  estimate 
the  complex  prediction  filter  coefficients  via  equation  (6)  which  was  solved 
using  an  efficient  form  of  the  Levinson  algorithm.  These  coefficients  were 
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augmented  by  128-Nj^  zeros  to  allow  the  FFT  routine  to  be  used  to  efficiently 

evaluate  the  denominator  of  equation  (5)  at  128  equispaced  points  in  wavenumber 
space  for  each  frequency.  The  number  of  wavenumber  bins  corresponding  to  plane 
waves  arriving  from  real  arrival  angles  varies  linearly  from  zero  at  d.c.  to 
128  at  fj^  (the  frequency  corresponding  to  the  half  wavelength  of  the  array). 

The  region  which  does  not  correspond  to  real  arrival  directions  is  termed  the 
'non-physical'  region  of  the  frequency  wavenumber  spectrum  and  a  full  discussion 
and  interpretation  of  this  region  is  given  at  reference  5. 

A 

The  power  spectrum  was  scaled  by  the  maximum  value  of  S|^(f,k)  and  displayed  in 
a  waterfall  format  in  figure  3  for  Nj^=16  and  figure  4  for  Nj^=24.  (The  32  hydro¬ 
phone  array  provides  a  maximum  of  32  lags) . 

For  comparison  the  Fourier  estimate  of  the  frequency  wavenumber  power  spectrum 
was  calculated  as  indicated  in  figure  5.  After  FFT'ing  the  acoustic  data  from 
each  hydrophone  to  the  same  frequency  resolution  as  for  the  MEM  a  second  FFT  was 
used  to  evaluate  the  spatial  summation  in  equation  (3)  (see  reference  5) .  Note 
that  since  only  32  hydrophones  were  analysed  only  32  independent  wavenumber^  bins 
are  possible.  At  each  frequency  the  spatial  array  of  hydrophone  data  (ie  Xj(f) 

for  j  =  0,1, ...31)  was  augmented  with  96  zeros  to  allow  the  FFT  routine  to  be 
used  to  efficiently  evaluate  the  spatial  Fourier  series.  This  also  allowed  a 
ready  comparison  with  the  maximum  entropy  estimates  since  the  same  number  of 
wavenumber  points  were  estimated  at  each  frequency  for  both  methods . 

Comparing  figures  3  and  5  the  increased  wavenumber  (and  hence  angular)  resolution 
of  the  MEM  is  apparent.  Comparing  figures  3  and  4  where  the  number  of  spatial 
lags  used  has  been  increased  from  16  to  24  respectively,  two  conclusions  emerge: 

(1)  The  resolution  has  been  improved  by  increasing  the  number  of  lags;  and 

(2)  At  some  frequencies  (eg  the  lowest  two  frequencies  of  figure  4)  and  the  MEM 
is  starting  to  become  unstable. 

If  the  number  of  lags  is  increased  further  this  instability  becomes  more  pronounced. 


6.  CONCLUSIONS 

The  maximum  entropy  technique  outlined  in  this  paper  is  a  computationally 
efficient  way  to  increase  the  spatial  resolution  of  estimates  of  the  frequency 
wavenumber  power  spectrum  of  data  from  a  linear  array  of  sensors. 

The  effect  of  bias  in  the  estimate  of  the  spatial  covariance  function  which  is 
caused  by  windowing  can  cause  misleading  results.  A  simple  technique  for 
removing  this  bias  has  been  shown  to  be  effective  when  applying  MEM  to  real  data. 
The  technique  uses  a  subset  of  the  available  lags  of  the  spatial  covariance 
function;  unfortunately  it  becomes  unstable  as  the  number  of  lags  approaches  the 
maximum  possible  (ie  the  number  of  hydrophones  in  the  array).  Further  work 
using  information  theoretical  ideas  for  determining  the  number  of  lags  to  be  used 
is  in  progress. 
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